***************************************************************************
* Project: Measuring Media Freedom: An Item Response Theory Analysis      *
* Authors: Jonathan A. Solis (jsolis@aiddata.wm.edu) & Philip D. Waggoner *
* Journal: British Journal of Political Science (2020)                    *
* Task: Replicate tables 4, 5, & 6; figure 4 in BJPS Appendix             *
* Date: 01/27/2020                                                        *
* Software: Stata 15.1                                                    *
***************************************************************************

* This .do file provides code to estimate the models found in table 4, 5, and 6 as well
* as create figure 4 in Solis and Waggoner's (2020) Appendix. This covers replications
* for Whitten-Woodring (2009), Stier (2015), and Schoonvelde (2014) respectively. We 
* dedicate a section to each table below (Section I includes code to replicate figure 4, 
* which we derive from models in table 4). Each section/table also has its own corresponding   
* dataset. In the replication materials, we provide the .dta file (Stata dataset) for 
* each dataset as well as an equivalent .csv file. 
*
* Below we list the three (3) datasets you will need to execute the code in this .do file:
*
*		(1) `bjps_msf-apndx_data_jww'
*		(2) `bjps_msf-apndx_data_stier'
*		(3) `bjps_msf-apdnx_data_schoonvelde'

				***************************************************************
				***************************************************************
				***************************************************************
                *****-----------------------------------------------------*****
				***** SECTION I. Replication from Whitten-Woodring (2009) *****
				*****-----------------------------------------------------*****
				***************************************************************
				***************************************************************
				***************************************************************

**********************************************************************************************
* Note: -We obtained the original dataset and Stata .do file of Whitten-Woodring's (2009)    *
* 		 "Watchdog or Lapdog? Media Freedom, Regime Type, and Government Respect for Human   *
*        Rights" and merged our new MSF point estimates and standard deviations              *
*       -Dataset obtained via Harvard's Dataverse:                                           *
*       --->https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/AGZHQW  *
**********************************************************************************************

*Load data
clear all
cd " "
use bjps_msf-apndx_data_jww
set more off

*Create interaction term with new MSF data
gen interaction= msf*FivePol

*Set data to time series
tsset ccode year, yearly

****************************************
*--------------------------------------*
******** TABLE 4, BJPS APPENDIX ********
*--------------------------------------*
****************************************

*Model 1, BJPS appendix (Table 4)
* -> original Whitten-Woodring model 8 in table 3 (pg. 612)
xtpcse physint Dmedfree FivePol FPolMF lnrgdp lnpop T2Conflict CConflict, correlation (ar1)

*Model 2, BJPS appendix (Table 4)
* -> original Whitten-Woodring model 8 in table 3 (pg. 612) but replace Freedom 
*    House data w/ MSF point estimates
xtpcse physint msf FivePol interaction lnrgdp lnpop T2Conflict CConflict, correlation (ar1)

**********************************************************
*Monte Carlo SIMULATION, BJPS APPENDIX (table 4, model 3)*
**********************************************************

****************************************************************************************
* Note: -Monte Carlo procedure requires us to generate new data (750 random draws      *
*        of MSF distribution) for each model                                           *
*       -Code repurposed from http://www.unified-democracy-scores.org/example.html     *
****************************************************************************************

*Model 3, BJPS appendix (Table 4)
* -> original Whitten-Woodring model 8 in table 3 (pg. 612) but replace Freedom 
*    House data w/ random draws of MSF posterior distribution
* -> runs model 750 times with values randomly pulled from the MSF posterior distribution

*Generate 750 random draws from MSF distribution *
* * * * * * * * * * * * * * * * * *  * * * * * * *

*Set seed for reproducibility
set seed 159484

*Loop for 750 draws
forvalues i = 1/750 {
gen x`i'=rnormal(msf,postsd)
}
*

*****************************
**** Prep for the Monte Carlo
*****************************

      set more off
      set matsize 750
	  
      *** Note that the matsize must be at least as large as the larger
      *** dimension in your posterior matrix. Since we use only 750 draws,
	  *** setting the matsize to 750 is sufficient for Stata IC, MP, and SE.
	  *** Executing a model with more draws where the dimensions rise above 800
	  *** will require a Stata more powerful than IC.

      *** Run the Monte Carlo
      forvalues i = 1/750 {
      *** Print out an iteration number
      display `i'

      *** Fit the model, using the ith draw from the MSF posterior
	  quietly xtpcse physint c.x`i'##o.FivePol lnrgdp lnpop T2Conflict CConflict, ///
      correlation (ar1)

	  *** Extract the coefficients and variance-covariance matrix
      matrix b = e(b)
      matrix V = e(V)
      local blength = colsof(b)

      *** Preserve the dataset, take a single multivariate normal draw from the
      *** posterior distribution of the coefficients, and restore the dataset.
      *** We use the capture command to catch possible errors in drawnorm
      *** and drop these iterations gracefully
      preserve
      capture quietly drawnorm b1-b`blength', double n(1) means(b) cov(V) clear
      if _rc == 0 {
        mkmat b1-b`blength', matrix(bsample)
        matrix posterior = nullmat(posterior) \ bsample
                  }
      else {
            display "Error drawing sample...iteration dropped"
           }
		   
      restore
	   } 
*

*** Create variables from matrix
svmat posterior

*** Calculate means and standard deviations
* -> NOTE: These are the results found in appendix table 4 (model 3)
tabstat posterior*, stat(mean sd)

*** Reset matrix size (not necessary to execute code but here if you need to reset the matrix size)
matrix drop _all

*****************************************
*---------------------------------------*
******** FIGURE 4, BJPS APPENDIX ********
*---------------------------------------*
*****************************************

*****************************************************************************			 
*Figure 4(a), BJPS APPENDIX                                                 *
* -> original Whitten-Woodring figure A4 (pg. 622)                           *
* -> code repurposed from: http://mattgolder.com/interactions               *
*****************************************************************************

*Run model and generate sample identifier for figure histogram
xtpcse physint Dmedfree FivePol FPolMF lnrgdp lnpop T2Conflict CConflict, correlation (ar1)
*Identify sample
gen byte used=e(sample)

*Pull coefficient vector
matrix b=e(b)
*Pull covariance matrix
matrix V=e(V)

*See coefficient vector
matrix list b
*See covariance matrix
matrix list V

*Save scalar of beta coefficient for "Dmedfree"
scalar b1=b[1,1]
*Save scalar of interaction "FPolMF"
scalar b3=b[1,3]

*Save scalar of variance of "Dmedfree" w/ itself
scalar varb1=V[1,1]
*Save scalar of variance of interaction "FPolMF" w/ itself
scalar varb3=V[3,3]

*Make scalar of covariance of X 'Dmedfree" and interaction "FPolMF"
scalar covb1b3=V[1,3]

*List of what you just created
scalar list b1 b3 varb1 varb3 covb1b3

*********************************************************
* Calculate data necessary for top marginal effect plot *
*********************************************************

*Generate variable MVZ that takes on all of the values of the modifying
*variable Z for which you want to calculate the marginal effect of X 
gen MVZ=((_n-.001)/100)

*Determine the last value of the modifying variable Z for which 
*you want to calculate the marginal effect of X. 
replace  MVZ=. if _n>401

*Calculate the marginal effect of X for all values MVZ,
*so long as MVZ is less than 402
gen conbx=b1+b3*MVZ if _n<402

*Calculate the standard error for the marginal effect of X 
*for all values MVZ, so long as MVZ is less than 402
gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ) if _n<402

*Use the standard errors and marginal effects to create two-tailed 95% confidence 
*intervals. These lines create the necessary information if you want to include
*a rug plot in the marginal effect plot (though we do not).
gen ax=1.96*consx
gen upperx=conbx+ax
gen lowerx=conbx-ax

******************************************
* Construct variable to produce y=0 line *
******************************************

*This will create a horizontal line at 0 on the vertical axis.
gen yline=0

**************************************
* Produce marginal effect plot for X *
**************************************
graph twoway hist FivePol if used==1, width(0.01) percent color(gs14) yaxis(2) ///
	    ||   line conbx   MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1) ///
        ||   line upperx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line lowerx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line yline  MVZ,  clwidth(thin) clcolor(black) clpattern(solid) ///
	    ||   , ///
             xlabel(0 1 2 3 4, nogrid labsize(2)) ///
		     ylabel(-2 -1 0 1 2, axis(1) nogrid labsize(2)) ///
		     ylabel(0 5 10 15 20 25 30 35 40, axis(2) nogrid labsize(2)) ///
	         yscale(noline alt) ///
		     yscale(noline alt axis(2)) ///	
             xscale(noline) ///
             legend(off) ///
			 subtitle("") ///
             xtitle("Regime Characteristics" , size(2.5)  ) ///
			 ytitle("Marginal Effect of Media Freedom" , axis(1) size(2.5)) ///
             ytitle("Percentage of Observations" , axis(2) size(2.5) orientation(rvertical)) ///
             xsca(titlegap(2)) ///
             ysca(titlegap(2)) ///
			 scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

*****************************************************************************			 
*Figure 4(b), BJPS APPENDIX                                                 *
* -> original Whitten-Woodring figure A4 (pg. 622) but replace Freedom House *
*    variable with MSF point estimates                                      *
* -> code repurposed from: http://mattgolder.com/interactions               *
*****************************************************************************

* <-> Drop variables generated from figure 4a or code won't run <->			 
drop used MVZ conbx consx ax upperx lowerx yline

*Run model and generate sample identifier for figure histogram
xtpcse physint msf FivePol interaction lnrgdp lnpop T2Conflict CConflict, correlation (ar1)
*Identify sample
gen byte used=e(sample)

*Pull coefficient vector
matrix b=e(b)
*Pull covariance matrix
matrix V=e(V)

*See coefficient vector
matrix list b
*See covariance matrix
matrix list V

*Save scalar of beta coefficient for "msf"
scalar b1=b[1,1]
*save scalar of interaction "interaction"
scalar b3=b[1,3]

*Save scalar of variance of "msf" w/ itself
scalar varb1=V[1,1]
*Save scalar of variance of interaction "interaction" w/ itself
scalar varb3=V[3,3]

*Make scalar of covariance of X 'msf" and interaction "interaction"
scalar covb1b3=V[1,3]

*List of what you just created
scalar list b1 b3 varb1 varb3 covb1b3

*********************************************************
* Calculate data necessary for top marginal effect plot *
*********************************************************

*Generate variable MVZ that takes on all of the values of the modifying
*variable Z for which you want to calculate the marginal effect of X 
gen MVZ=((_n-.001)/100)

*Determine the last value of the modifying variable Z for which 
*you want to calculate the marginal effect of X. 
replace  MVZ=. if _n>401

*Calculate the marginal effect of X for all values MVZ,
*so long as MVZ is less than 402
gen conbx=b1+b3*MVZ if _n<402

*Calculate the standard error for the marginal effect of X 
*for all values MVZ, so long as MVZ is less than 402
gen consx=sqrt(varb1+varb3*(MVZ^2)+2*covb1b3*MVZ) if _n<402

*Use the standard errors and marginal effects to create two-tailed 95% confidence 
*intervals. These lines create the necessary information if you want to include
*a rug plot in the marginal effect plot (though we do not).
gen ax=1.96*consx
gen upperx=conbx+ax
gen lowerx=conbx-ax

******************************************
* Construct variable to produce y=0 line *
******************************************

*This will create a horizontal line at 0 on the vertical axis.
gen yline=0

**************************************
* Produce marginal effect plot for X *
**************************************
graph twoway hist FivePol if used==1, width(0.01) percent color(gs14) yaxis(2) ///
	    ||   line conbx   MVZ, clpattern(solid) clwidth(medium) clcolor(black) yaxis(1) ///
        ||   line upperx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line lowerx  MVZ, clpattern(dash) clwidth(thin) clcolor(black) ///
        ||   line yline  MVZ,  clwidth(thin) clcolor(black) clpattern(solid) ///
	    ||   , ///
             xlabel(0 1 2 3 4, nogrid labsize(2)) ///
		     ylabel(-4 -2 0 2 4 6 8, axis(1) nogrid labsize(2)) ///
		     ylabel(0 5 10 15 20 25 30 35 40, axis(2) nogrid labsize(2)) ///
	         yscale(noline alt) ///
		     yscale(noline alt axis(2)) ///	
             xscale(noline) ///
             legend(off) ///
			 subtitle("") ///
             xtitle("Regime Characteristics" , size(2.5)  ) ///
			 ytitle("Marginal Effect of Media Freedom" , axis(1) size(2.5)) ///
             ytitle("Percentage of Observations" , axis(2) size(2.5) orientation(rvertical)) ///
             xsca(titlegap(2)) ///
             ysca(titlegap(2)) ///
			 scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white))

*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*

				*****************************************************
				*****************************************************
				*****************************************************
                *****-------------------------------------------*****
				***** Section II. Replication from Stier (2015) *****
				*****-------------------------------------------*****
				*****************************************************
				*****************************************************
				*****************************************************

****************************************************************************************
* Note: -We obtained the original dataset and Stata code of Stier et al.'s (2015)      *
* 		 "Democracy, Autocracy and the News: The Impact of Regime Type on Media        *
*        Freedom" and merged our new MSF point estimates and standard deviations       *
*       -Dataset obtained via email from article author Sebastian Stier Nov 5, 2018    *
****************************************************************************************

*Load data
clear all
cd " "
use bjps_msf-apndx_data_stier
set more off

*Clean
drop if cowcode==.

*Set to time-series
tsset cowcode year

****************************************
*--------------------------------------*
******** TABLE 5, BJPS APPENDIX ********
*--------------------------------------*
****************************************

*Model 4, BJPS appendix (Table 5) 
* -> original Stier model 1 in table 1 (pg. 1,284)
* -> syntax "gen byte m1_1=e(sample)" creates variable to Identify sample
*    used in this model; we use this to create a comparable sample with
*    the new MSF data model
xtpcse fopneu l.kai_libdem l.log_gdp l.log_pop l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_transition==0 ///
 & kai_fail ==0, correlation(psar1) pairwise
*Identify sample
gen byte m1_1=e(sample)

*Model 5, BJPS appendix (Table 5) 
* -> original Stier model 1 in table 1 (pg. 1,284) but replace Freedom 
*    House data w/ MSF point estimates
* -> syntax "m1_1==1" ensures we analyze a comparable sample w/ new MSF data
xtpcse msf l.kai_libdem l.log_gdp l.log_pop l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_transition==0 ///
 & kai_fail ==0 & m1_1==1, correlation(psar1) pairwise

*Model 6, BJPS appendix (Table 5) 
* -> see Monte Carlo Simulations section below

*Model 7, BJPS appendix (Table 5) 
* -> original Stier model 3 in table 1 (pg. 1,284)
* -> syntax "gen byte m3_1=e(sample)" creates variable to Identify sample
*    used in this model; we use this to create a comparable sample with
*    the new MSF data model
xtpcse fopneu l.kai_monarchy l.kai_military l.kai_oneparty l.kai_personalist ///
 l.kai_communist l.log_gdp l.log_pop  l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_libdem==0 ///
 & kai_transition==0 & kai_fail ==0 & cowcode != 553 & cowcode != 800, ///
 correlation(psar1) pairwise
*Identify sample
gen byte m3_1=e(sample)

*Model 8, BJPS appendix (Table 5) 
* -> original Stier model 3 in table 1 (pg. 1,284) but replace Freedom 
*    House data w/ MSF point estimates
* -> syntax "m3_1==1" ensures we analyze a comparable sample w/ new MSF data
xtpcse msf l.kai_monarchy l.kai_military l.kai_oneparty l.kai_personalist ///
 l.kai_communist l.log_gdp l.log_pop  l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_libdem==0 ///
 & kai_transition==0 & kai_fail ==0 & cowcode != 553 & cowcode != 800  ///
 & m3_1==1,  correlation(psar1) pairwise

*Model 9, BJPS appendix (Table 5) 
* -> see Monte Carlo Simulations section below

*****************************************************************
* Monte Carlo SIMULATIONS, BJPS Appendix Table 5 (models 6 & 9) *
*****************************************************************

****************************************************************************************
* Note: -To execute the Monte Carlo Simulations, we clear and reload original dataset  *
*       -Monte Carlo procedure requires us to generate new data (750 random draws      *
*        of MSF lead distribution) for each model                                      *
*       -Code repurposed from http://www.unified-democracy-scores.org/example.html     *
****************************************************************************************

*Model 6, BJPS appendix (Table 5) 
* -> original Stier, model 1 in table 1 (pg. 1,284) but replace Freedom 
*    House data w/ MSF point estimates
* -> runs model 750 times with values randomly pulled from the MSF distribution
*    in each of the 750 models 

*Load data
clear all
cd " "
use bjps_msf-apndx_data_stier
set more off

*Clean
drop if cowcode==.

*Set to time series
tsset cowcode year

*Run model to mark observations in Stier's original sample for comparability
quietly xtpcse fopneu l.kai_libdem l.log_gdp l.log_pop l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_transition==0 ///
 & kai_fail ==0, correlation(psar1) pairwise
*Identify sample
 gen byte m1_1=e(sample)

*Generate 750 random draws from MSF distribution *
* * * * * * * * * * * * * * * * * * * * ** * * * *

*Set seed for reproducibility
set seed 58684828

*Loop for 750 draws
forvalues i = 1/750 {
gen x`i'=rnormal(msf,postsd)
}
*

*****************************
**** Prep for the Monte Carlo
*****************************

      set more off
      set matsize 750
	  
      *** Note that the matsize must be at least as large as the larger
      *** dimension in your posterior matrix. Since we use only 750 draws,
	  *** setting the matsize to 750 is sufficient for Stata IC, MP, and SE.
	  *** Executing a model with more draws where the dimensions rise above 800
	  *** will require a Stata more powerful than IC.
	  
      *** Run the Monte Carlo
      forvalues i = 1/750 {
      *** Print out an iteration number
      display `i'

      *** Fit the model, using the ith draw from the MSF posterior
	  quietly xtpcse x`i' l.kai_libdem l.log_gdp l.log_pop l.ross_oil_gas_valuepop_2009 ///
      l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_transition==0 ///
      & kai_fail ==0 & m1_1==1, correlation(psar1) pairwise

	  *** Extract the coefficients and variance-covariance matrix
      matrix b = e(b)
      matrix V = e(V)
      local blength = colsof(b)

      *** Preserve the dataset, take a single multivariate normal draw from the
      *** posterior distribution of the coefficients, and restore the dataset.
      *** We use the capture command to catch possible errors in drawnorm
      *** and drop these iterations gracefully.
      preserve
      capture quietly drawnorm b1-b`blength', double n(1) means(b) cov(V) clear
      if _rc == 0 {
        mkmat b1-b`blength', matrix(bsample)
        matrix posterior = nullmat(posterior) \ bsample
                  }
      else {
           display "Error drawing sample...iteration dropped"
           }
		   
      restore
      
	   } 
*

*** Create variables from matrix
svmat posterior

*** Calculate means and standard deviations
* -> NOTE: These are the results found in table 5 (model 6, BJPS appendix)
tabstat posterior*, stat(mean sd)

*** Reset matrix size (not necessary to execute code but here if you need to reset the matrix size)
matrix drop _all

*Model 9, BJPS appendix (Table 5)
* -> original Stier, model 3 in table 1 (pg. 1,284) but replace Freedom 
*    House data w/ MSF point estimates
* -> runs model 750 times with values randomly pulled from the MSF distribution
*    in each of the 750 models 

*Load data
clear all
cd " "
use bjps_msf-apndx_data_stier
set more off

*Clean
drop if cowcode==.

*Set to time-series
tsset cowcode year

*Run model to mark observations in Stier's original sample for comparability
quietly xtpcse fopneu l.kai_monarchy l.kai_military l.kai_oneparty l.kai_personalist ///
 l.kai_communist l.log_gdp l.log_pop  l.ross_oil_gas_valuepop_2009 ///
 l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_libdem==0 ///
 & kai_transition==0 & kai_fail ==0 & cowcode != 553 & cowcode != 800, ///
 correlation(psar1) pairwise
*Identify sample
 gen byte m3_1=e(sample)

*Generate 750 random draws from MSF distribution *
* * * * * * * * * * * * * * * * * * * * ** * * * *

*Set seed for reproducibility
set seed 19824

*Loop for 750 draws
forvalues i = 1/750 {
gen x`i'=rnormal(msf,postsd)
}
*

*****************************
**** Prep for the Monte Carlo
*****************************

      set more off
      set matsize 750
	  
      *** Note that the matsize must be at least as large as the larger
      *** dimension in your posterior matrix. Since we use only 750 draws,
	  *** setting the matsize to 750 is sufficient for Stata IC, MP, and SE.
	  *** Executing a model with more draws where the dimensions rise above 800
	  *** will require a Stata more powerful than IC.
	  
      *** Run the Monte Carlo
      forvalues i = 1/750 {
      *** Print out an iteration number
      display `i'

      *** Fit the model, using the ith draw from the MSF posterior
      quietly xtpcse x`i' l.kai_monarchy l.kai_military l.kai_oneparty l.kai_personalist ///
      l.kai_communist l.log_gdp l.log_pop  l.ross_oil_gas_valuepop_2009 ///
      l.inetusers_ipo l.incidencev412 i.own_aclp d00 if kai_libdem==0 ///
      & kai_transition==0 & kai_fail ==0 & cowcode != 553 & cowcode != 800 ///
      & m3_1==1, correlation(psar1) pairwise


	  *** Extract the coefficients and variance-covariance matrix
      matrix b = e(b)
      matrix V = e(V)
      local blength = colsof(b)

      *** Preserve the dataset, take a single multivariate normal draw from the
      *** posterior distribution of the coefficients, and restore the dataset.
      *** We use the capture command to catch possible errors in drawnorm
      *** and drop these iterations gracefully.
      preserve
      capture quietly drawnorm b1-b`blength', double n(1) means(b) cov(V) clear
      if _rc == 0 {
        mkmat b1-b`blength', matrix(bsample)
        matrix posterior = nullmat(posterior) \ bsample
                  }
      else {
           display "Error drawing sample...iteration dropped"
           }
      restore
      
	   } 
*

***** Create variables from matrix
svmat posterior

*** Calculate means and standard deviations
* -> NOTE: These are the results found in table 5 (model 9, BJPS appendix)
tabstat posterior*, stat(mean sd)

*** Reset matrix size (not necessary to execute code but here if you need to reset the matrix size)
matrix drop _all

*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*
*---------------------------------------------------------------------------------------------------*

				************************************************************
				************************************************************
				************************************************************
                *****--------------------------------------------------*****
				***** Section III. Replication from Schoonvelde (2014) *****
				*****--------------------------------------------------*****
				************************************************************
				************************************************************
				************************************************************

**********************************************************************************************
* Note: -We obtained the original dataset of Schoonvelde's (2014) "Media Freedom and the     *
* 		 Institutional Underpinnings of Political Knowledge" and merged our new MSF point    *
*        estimates and standard deviations                                                   *
*       -Dataset obtained via Harvard's Dataverse:                                           *
*       --->https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/24122   *
*       -The author performed the original analysis in R; However, we replicate the          *
*        the regression models in Stata below in order to perform the Monte Carlo simulation *
*        code which was written in Stata                                                     *    
**********************************************************************************************

*Load data
clear all
cd " "
use bjps_msf-apdnx_data_schoonvelde
set more off

*Set data to panel (cyear is the combined COW code and year; e.g., Australia's 
*cyear value for election year 1996 is '9001996'
*NOTE: Belgium 
xtset cyear

****************************************
*--------------------------------------*
******** TABLE 6, BJPS APPENDIX ********
*--------------------------------------*
****************************************

*Model 10, BJPS appendix (Table 6)
* -> original Schoonvelde model 2 in table 4 (pg. 172)
xtreg sophistication income01_wi_01 education01_wi_01 union age ///
freedom01 parliamentary compvoting partylist majoritarian enp herfindahln ///
herfindahlt newspaper public gdp

*Model 11, BJPS appendix (Table 6)
* -> original Schoonvelde model 2 in table 4 (pg. 172) but replace Freedom 
*    House data w/ MSF point estimates
xtreg sophistication income01_wi_01 education01_wi_01 union age ///
msf parliamentary compvoting partylist majoritarian enp herfindahln ///
herfindahlt newspaper public gdp

*Model 12, BJPS appendix (Table 6) 
* -> see Monte Carlo Simulations section below

*Model 13, BJPS appendix (Table 6)
* -> original Schoonvelde model 4 in table 4 (pg. 172)
xtreg sophistication income01_wi_01 education01_wi_01 freedom01 ///
parliamentary compvoting partylist majoritarian gdp

*Model 14, BJPS appendix (Table 6)
* -> original Schoonvelde model 4 in table 4 (pg. 172) but replace Freedom 
*    House data w/ MSF point estimates
xtreg sophistication income01_wi_01 education01_wi_01 msf ///
parliamentary compvoting partylist majoritarian gdp

*Model 15, BJPS appendix (Table 6) 
* -> see Monte Carlo Simulations section below

*******************************************************************
* Monte Carlo SIMULATIONS, BJPS Appendix Table 6 (models 12 & 15) *
*******************************************************************

****************************************************************************************
* Note: -To execute the Monte Carlo Simulations, we clear and reload original dataset  *
*       -Monte Carlo procedure requires us to generate new data (750 random draws      *
*        of MSF lead distribution) for each model                                      *
*       -Code repurposed from http://www.unified-democracy-scores.org/example.html     *
****************************************************************************************

*Model 12, BJPS appendix (Table 6)
* -> original Schoonvelde, model 2 in table 4 (pg. 172) but replace Freedom House 
*    data w/ MSF point estimates
* -> runs model 750 times with values randomly pulled from the MSF distribution
*    in each of the 750 models 

*Load data
clear all
cd " "
use bjps_msf-apdnx_data_schoonvelde
set more off

*Set to panel
xtset cyear

*Generate 750 random draws from MSF distribution*
* * * * * * * * * * * * * * * * * * * * * * * * *

*Set seed for reproducibility
set seed 159617

*Loop for 750 draws
forvalues i = 1/750 {
gen x`i'=rnormal(msf,postsd)
}
*

*****************************
**** Prep for the Monte Carlo
*****************************

      set more off
      set matsize 750
	  
      *** Note that the matsize must be at least as large as the larger
      *** dimension in your posterior matrix. Since we use only 750 draws,
	  *** setting the matsize to 750 is sufficient for Stata IC, MP, and SE.
	  *** Executing a model with more draws where the dimensions rise above 800
	  *** will require a Stata more powerful than IC.
	  
      *** Run the Monte Carlo
      forvalues i = 1/750 {
        *** Print out an iteration number
        display `i'

      *** Fit the model, using the ith draw from the MSF posterior	
	  quietly xtreg sophistication income01_wi_01 education01_wi_01 union age ///
	  x`i' parliamentary compvoting partylist majoritarian enp herfindahln ///
	  herfindahlt newspaper public gdp

	   *** Extract the coefficients and variance-covariance matrix
       matrix b = e(b)
       matrix V = e(V)
       local blength = colsof(b)

       *** Preserve the dataset, take a single multivariate normal draw from the
       *** posterior distribution of the coefficients, and restore the dataset.
       *** We use the capture command to catch possible errors in drawnorm
       *** and drop these iterations gracefully.
       preserve
       capture quietly drawnorm b1-b`blength', double n(1) means(b) cov(V) clear
       if _rc == 0 {
         mkmat b1-b`blength', matrix(bsample)
         matrix posterior = nullmat(posterior) \ bsample
                   }
        else {
          display "Error drawing sample...iteration dropped"
             }
        restore
      
	          }  
*

*** Create variables from matrix
svmat posterior

*** Calculate means and standard deviations
* -> NOTE: These are the results found in table 6 (model 12, BJPS appendix)
tabstat posterior*, stat(mean sd)

*** Reset matrix size (not necessary to execute code but here if you need to reset the matrix size)
matrix drop _all

*Model 15, BJPS appendix (Table 6)
* -> original Schoonvelde, model 4 in table 4 (pg. 172) but replace Freedom House 
*    data w/ MSF point estimates
* -> runs model 750 times with values randomly pulled from the MSF distribution
*    in each of the 750 models 

*Load data
clear all
cd " "
use bjps_msf-apdnx_data_schoonvelde
set more off

*Set to panel
xtset cyear

*Generate 750 random draws from MSF distribution *
* * * * * * * * * * * * * * * * * ** * * * * * * *

*Set seed for reproducibility
set seed 716951

*Loop for 750 draws
forvalues i = 1/750 {
gen x`i'=rnormal(msf,postsd)
}
*

*****************************
**** Prep for the Monte Carlo
*****************************

      set more off
      set matsize 750
	  
      *** Note that the matsize must be at least as large as the larger
      *** dimension in your posterior matrix. Since we use only 750 draws,
	  *** setting the matsize to 750 is sufficient for Stata IC, MP, and SE.
	  *** Executing a model with more draws where the dimensions rise above 800
	  *** will require a Stata more powerful than IC.
	  
      *** Run the Monte Carlo
      forvalues i = 1/750 {
      *** Print out an iteration number
      display `i'
	
	  *** Fit the model, using the ith draw from the MSF posterior
	  quietly xtreg sophistication income01_wi_01 education01_wi_01 x`i' ///
	  parliamentary compvoting partylist majoritarian gdp
 
      *** Extract the coefficients and variance-covariance matrix
      matrix b = e(b)
      matrix V = e(V)
      local blength = colsof(b)

      *** Preserve the dataset, take a single multivariate normal draw from the
      *** posterior distribution of the coefficients, and restore the dataset.
      *** We use the capture command to catch possible errors in drawnorm
      *** and drop these iterations gracefully.
      preserve
      capture quietly drawnorm b1-b`blength', double n(1) means(b) cov(V) clear
      if _rc == 0 {
         mkmat b1-b`blength', matrix(bsample)
         matrix posterior = nullmat(posterior) \ bsample
                   }
      else {
          display "Error drawing sample...iteration dropped"
            }
      restore
      
	     }  
*

*** Create variables from matrix
svmat posterior

*** Calculate means and standard deviations
* -> NOTE: These are the results found in table 6 (model 15, BJPS appendix)
tabstat posterior*, stat(mean sd)

*** Reset matrix size (not necessary to execute code but here if you need to reset the matrix size)
matrix drop _all

*-----------------------------------------------------------------------------*
*-----------------------------------------------------------------------------*
*----------------------------------*********----------------------------------*
*----------------------------------** END **----------------------------------*
*----------------------------------*********----------------------------------*
*-----------------------------------------------------------------------------*
*-----------------------------------------------------------------------------*
